(1) Field of the Invention
This invention relates to time-domain simulation of the wave equation and, in particular, to a method and a system to compute an evolution of the wave equation.
(2) Description of Related Art
In one dimension, the derivative of a point on a line is the slope of that line at that point. The slope at any point on a line is the tangent line through a particular point, which has the same direction as the line at that point. In other words, the derivative of a point on a line is construed as a local statement about the local change in the direction of the line. The partial differential equation of a function with four variables is defined in much the same way, by using a tangent to a three dimensional function in a way analogous to the way tangent line is used for a single variable function. Partial differential equations may be used to represent a particular point of a field in space and time. The slope (derivative) of a field at which any one point exists depends on the three variables for space in x, y, and z directions and one variable for time t. Therefore, partial differential equations can be construed as local statements about local changes in a three dimensional entity, such as a wave, in space and time.
The wave Equation (1) is represented by a partial differential equation because it involves a function of four variables, which are the three dimensions of space in the x, y, and z directions and one variable for time t as follows,
                                                        (                                                ∇                  2                                ⁢                                  -                                                            ∂                      2                                                              ∂                                              t                        2                                                                                                        )                        ⁢                          ψ              ⁡                              (                                                      x                    ρ                                    ,                  t                                )                                              =          0                ,                            (        1        )            where:
      ∇    2    ⁢      =                            ∂          2                          ∂                      x            2                              +                        ∂          2                          ∂                      y            2                              +                        ∂          2                          ∂                      z            2                              denotes the well-known Laplace operator,ψ(, t) is the field quantity, and=dimensions of space in x, y, and z directions.
Wave fields are continuous entities in both space and time, which means that they have a value at every point in space and time. This makes the use of partial differential equations for an exact representation of a field at every point in space and time impractical and difficult because an infinite number of points must be used to provide an exact representation of the total field at every point in space and time. Other analyzing techniques must therefore be used to avoid the use of partial differential equations for representing a continuous entity such as a wave.
Numerical Analysis is one method for computing a solution for partial differential equations. This is done by first discretizing a partial differential equation, bringing it into a finite dimensional subspace, then solving the resulting linear system in the finite dimensional subspace. One such technique for first discretizing the wave Equation (1) involves using a finite differencing method. This method approximates the solution of partial differential equations by using a finite collection of points to represent a wave. The approach taken for partial differential equations is to approximate differential operators such as u′(x) by a difference operator such as (u(x+h)−u(x))/h for some small but finite h, where h may represent a small rate of change Δ. Doing this substitution for a sufficient number of points in the domain of definition (for instance 0, h, 2h . . . 1 in the case of the unit interval) provides a system of equations that can be solved algebraically.
The process of representing a continuous entity by a finite collection of points is known as discretization. In other words, discretization is a way of approximating a problem that has an infinite number of variables with a problem that has a finite number of variables. A discretization process divides geometry into smaller pieces (finite elements) to prepare for analysis. A discretized model of a continuous entity such as a wave can be manipulated using an analyzing technique such as a Finite-Difference Time-Domain (FDTD), in which all the simple geometric components of a model are evaluated simultaneously. To apply the FDTD, the differential equation of the wave is replaced by a finite difference equation. This is done by replacing the partial derivatives by central difference quotients, which replaces partial differential equations with algebraic expressions. In the next operations, the interval to be computed is subdivided into a suitable number of equal subintervals, and then the difference equation at each point where the value of the function is unknown is determined. At the final operation, the system of equations created is solved to obtain approximation values for the solution of the differential equation at discrete points on the interval.
Using FDTD, at a given position in space, the field at each current step can be calculated from the previous values of the fields (hence the difference in time). If time is considered discretized into discrete steps, then fields are simply increased or decreased from their previously calculated values by adding or subtracting an amount given by the
calculations at current time. The partial differential equations, which describe the infinitesimal local behavior of the field, are replaced with algebraic equations having to do with differences in the field in nearby space and time points. This avoids the partial differential operations because the equations are transformed into algebraic expressions. Dividing a real problem into appropriately-sized unit cells, and then calculating the field equations, iterated for as long as the response is of interest, therefore solves a real problem.
To use FDTD, a computational domain, which is the space in which the simulation of a wave will be performed, must be established. The field will be determined at every point within the computational domain, requiring a cell or grid structuring of the domain. In addition, since a computational domain must end at some point, a boundary must be established. Once the computational domain is determined, a source is then specified. The source can be an impinging plane wave or a current on a wire, depending on the type of situation to be modeled.
Since FDTD requires that the entire computational domain be formed of grids, and because these grids must be small compared to the smallest wavelength and smaller than the smallest feature model, very large computational domains can be developed, which result in very long solution times. Models with long, thin features are difficult to model in FDTD because of the excessively large computational domain required. In addition, although fields can be found directly everywhere in the computational domain when using FDTD, if, however, calculations or models of fields at some distance are desired, it is likely that the distance will force the computational domain to be excessively large. It should be noted that far-field extensions are available for FDTD, but require some amount of post-processing.
As is described in the Published Patent Application 2002/0099510 to Namiki, portions of which are repeated herein for completeness, conventionally, FDTD solutions are numerically obtained through explicit methods. For simplicity, consider the following one-dimensional parabolic partial differential equation,
                                          ∂            f                                ∂            t                          =                                                            ∂                2                            ⁢              f                                      ∂                              x                2                                              .                                    (        2        )            Equation (2) may be discretized with central difference approximations taken for the derivatives of time t and space x,ƒ(t,x)=ƒ(nΔt, iΔx)=ƒin.  (3)The Equation (3) indicates that ƒ is a function of time t and position x, where t and x are discretized as nΔt and iΔx, respectively. The temporal steps are indexed by some integer n and related to the continuous time by the relation t=nΔt, and spatial steps are indexed by another integer i and related to continuous space (e.g., in the x direction) by the relation x=iΔx. Therefore, n is the number of time steps that have elapsed from the beginning of computation, Δt is the temporal discretization interval (i.e., time step size), i is the grid coordinate number representing a specific point in the one-dimensional space, and Δx is the spatial discretization interval (i.e., space increment or cell size). To yield a solution by using an explicit method, the above differential equation is approximated to the following forward-difference expression with respect to its time derivative,
                                                        f              i                              n                +                1                                      -                          f              i              n                                            Δ            ⁢                                                  ⁢            t                          =                                                            f                                  i                  +                  1                                n                            -                              2                ⁢                                  f                  i                  n                                            +                              f                                  i                  -                  1                                n                                                    Δ              ⁢                                                          ⁢                              x                2                                              .                                    (        4        )            Let
      r    =                  Δ        ⁢                                  ⁢        t                    Δ        ⁢                                  ⁢                  x          2                      ,then Equation (4) can be rewritten as follows,ƒin+1=rƒi+1n+(1−2r)ƒin+rƒn−1n.  (5)Note that, when a solution at time nΔt is given, Equation (5) immediately gives the next solution at time (n+1)Δt. Numerical solvers of this kind are referred to as explicit methods.
The explicit difference method of the Equation (5), however, must satisfy the following condition for stability to make sure that the solution will converge toward a final solution,
                    r        =                                            Δ              ⁢                                                          ⁢              t                                      Δ              ⁢                                                          ⁢                              x                2                                              ≤                                    1              2                        .                                              (        6        )            Furthermore, to avoid numerical instability in FDTD computation, the following condition should be satisfied, which is known as the Courant, Friedrich, and Levy (CFL) condition, or Courant condition,
                              Δ          ⁢                                          ⁢          t                ≤                              1                                                                                                      (                                              1                                                  Δ                          ⁢                                                                                                          ⁢                                                      x                            min                                                                                              )                                        2                                    +                                                            (                                              1                                                  Δ                          ⁢                                                                                                          ⁢                                                      y                            min                                                                                              )                                        2                                                  ν                            +                                                (                                      1                                          Δ                      ⁢                                                                                          ⁢                                              z                        min                                                                              )                                2                                              .                                    (        7        )            
Detailed discussion on the CFL condition in FDTD is found in the literature, as described by A. Taflove, in “Computational Electrodynamics,” MA, Artech House Inc., 1995. The Equation (7) is known as the Courant condition for three-dimensional wave analysis, where v is the propagation rate of electromagnetic waves, and Δxmin, Δymin, and Δzmin are minimum values of spatial discretization intervals in the x, y, and z directions, respectively. Because of the constraints discussed above, the maximum time-step size in explicit methods is dependent on the minimum cell size. That is, the time-step size must be reduced when analyzing an object having fine geometrical features, resulting in an increased number of simulation steps to be iterated. This would cause a serious problem of long simulation time, particularly when calculating a time response for an extended period.
Partial differential equations can be solved with explicit as well as implicit methods. For example, a backward difference approximation to the one-dimensional parabolic partial differential Equation (2) is as follows,
                                                        f              i                              n                -                1                                      -                          f              i              n                                            Δ            ⁢                                                  ⁢            t                          =                                                            f                                  i                  +                  1                                                  n                  +                  1                                            -                              2                ⁢                                  f                  i                                      n                    +                    1                                                              +                              f                                  i                  -                  1                                                  n                  +                  1                                                                    Δ              ⁢                                                          ⁢                              x                2                                              .                                    (        8        )            
Similar to the case of explicit methods mentioned above, by allowing
      r    =                  Δ        ⁢                                  ⁢        t                    Δ        ⁢                                  ⁢                  x          2                      ,the difference Equation (8) can be rearranged as follows,−rƒi+1n+1+(1−2r)ƒin+1−rƒi−1n+1=ƒin.  (9)Equation (9) is an implicit expression of the problem to be solved. Unlike the explicit methods, the numerical stability is guaranteed when solving this implicit expression. In the implicit method, however, it is necessary to solve the following set of simultaneous equations in order to obtain a series of ƒin+1,
                                                        (                                                                                          b                      1                                                                                                  c                      1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            a                      2                                                                                                  b                      2                                                                                                  c                      2                                                                                                                                                                              0                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                ⋰                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      ⋰                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  0                                                                                                                                                                                a                                              imax                        -                        1                                                                                                                        b                                              imax                        -                        1                                                                                                                        c                                              imax                        -                        1                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  a                      imax                                                                                                  b                      imax                                                                                  )                        ⁢                          (                                                                                          f                      1                                              n                        +                        1                                                                                                                                                        f                      2                                              n                        +                        1                                                                                                                                  ⋮                                                                                                              f                                              imax                        -                        1                                                                    n                        +                        1                                                                                                                                                        f                      imax                                              n                        +                        1                                                                                                        )                                =                      (                                                                                f                    1                    n                                                                                                                    f                    2                    n                                                                                                ⋮                                                                                                  f                                          imax                      -                      1                                        n                                                                                                                    f                    imax                    n                                                                        )                          ,                            (        10        )            where ai is the invariable part (−r) of the third term of the left-hand side of Equation (9) when i=1, 2 . . . imax, bi is the invariable part (1+2r) of the second term of the left-hand side of Equation (9) when i=1, 2, . . . imax, and ci is the invariable part (−r) of the first term of the left-hand side of Equation (9) when i=1, 2, . . . imax.
To analyze the transient behavior of electromagnetic waves, it is required to solve a partial differential equation in the space of at least two dimensions. Special care must be taken when solving this kind of problem by use of implicit methods. Consider the following two-dimensional partial differential equation,
                                          ∂            f                                ∂            t                          =                                                            ∂                2                            ⁢              f                                      ∂                              x                2                                              +                                                                      ∂                  2                                ⁢                f                                            ∂                                  y                  2                                                      .                                              (        11        )            To solve Equation (11) with, for example, the well-known Crank-Nicolson method, its time-derivative term is approximated as follows:
                                                        f                              i                ,                j                                            n                +                1                                      -                          f                              i                ,                j                            n                                            Δ            ⁢                                                  ⁢            t                          =                                            1              2                        ⁢                                                            f                                                            i                      +                      1                                        ,                    j                                    n                                -                                  2                  ⁢                                      f                                          i                      ,                      j                                        n                                                  +                                  f                                                            i                      -                      1                                        ,                    j                                    n                                                            Δ                ⁢                                                                  ⁢                                  x                  2                                                              +                                                    f                                  i                  ,                                      j                    +                    1                                                  n                            -                              2                ⁢                                  f                                      i                    ,                    j                                    n                                            +                              f                                  i                  ,                                      j                    -                    1                                                  n                                                    Δ              ⁢                                                          ⁢                              y                2                                              +                                    1              2                        ⁢                                                            f                                                            i                      +                      1                                        ,                    j                                                        n                    +                    1                                                  -                                  2                  ⁢                                      f                                          i                      ,                      j                                                              n                      +                      1                                                                      +                                  f                                                            i                      -                      1                                        ,                    j                                                        n                    +                    1                                                                              Δ                ⁢                                                                  ⁢                                  x                  2                                                              +                                                                      f                                      i                    ,                                          j                      +                      1                                                                            n                    +                    1                                                  -                                  2                  ⁢                                      f                                          i                      ,                      j                                                              n                      +                      1                                                                      +                                  f                                      i                    ,                                          j                      -                      1                                                                            n                    +                    1                                                                              Δ                ⁢                                                                  ⁢                                  y                  2                                                      .                                              (        12        )            
The solution of Equation (12) can be reached by solving a set of simultaneous linear equations having as many unknowns as (Nx−1)×(Ny−1), where Nx and Ny are the numbers of meshes in the x-axis and y-axis directions, respectively. This computation makes extreme demands on the computer memory and processing power, even when fine meshing is required. The result is that implicit methods are as time-consuming as explicit methods, even when the simulation model has a fine geometrical feature that needs a smaller cell size.
There are two main issues with respect to performance of a numerical solution of the wave equations, stability and efficiency. The numerical stability of a method describes how that method responds to the differences between the calculation and the function being approximated. In a stable method, the errors due to the approximations are damped out as the computation proceeds. In an unstable method, any errors in processing are magnified as the calculation proceeds. Unstable methods quickly generate worthless data, and are useless for numerical processing. Therefore, any method used to compute the approximate solution of a wave equation must be stable.
The second issue is with respect to the efficiency in the numerical solutions of partial differential equations, which relates to the precision of the approximated solution. Recall that numerical solutions are mere approximations of the partial differential equations, and hence, the error between the approximate solution and the true solution must somehow be accounted for to determine the precision of the approximation. For example, is the approximation error small enough that it provides an acceptable solution for which the partial equations were modeled? In general, the error between the approximate solution and the true solution is determined by the truncation error that is made by going from a differential operator to a difference operator. The truncation error reflects the fact that a difference operator can be viewed as a finite part of the infinite series of the differential operator. Therefore, in general, a higher number of series (an increased number of steps or values closer together in space and time) can be used in the numerical solution to better approximate the partial differential equation to reduce the approximation error. In fact, the approximation error should reduce to zero if enough points in space and time are used to approximate the solution to the partial differential equations, which is the definition of convergence for an infinite series.
Assuming that the method is convergent, then a determination must be made as to how rapidly the approximation error reduces to zero (or the solution converges), which is the measure of the efficiency of the method being used. Most of the methods that are used today are low-order discretization methods, which means that the approximation error is decreased by a very low factor despite an increased use of points (numbers) in space and time. Discretizations with approximation errors that reduce rapidly as more points are used are known as high-order discretization methods.
The difficulty in changing a low-order time-domain simulation to a high-order one is that, in general, the low-order time-domain simulations tend to be unstable high-order solutions. Another difficulty is that in order to achieve a high-order solution, the algebraic equations that replace the partial differential equations would have to relate to the values of the field at many steps (time and or space). One major problem with this is initialization. For example, if the value of a previous four time steps are needed, then initial conditions for the four time steps are required, and not just one. On the other hand, the more time steps that are required, the greater the possibility of introducing extraneous and unstable solutions.
In light of the current state of the art and the drawbacks to current systems mentioned above, a need exists for a system and a method that would allow for stable, high-order discretization for evolution of the wave equation.